#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>
#include "matrix.h"

gsl_matrix* invertMatrix(gsl_matrix* matrix) {
	int s;
	int n = matrix->size1;
	gsl_matrix * inverse = gsl_matrix_alloc (n, n);
	gsl_permutation * perm = gsl_permutation_alloc (n);
	gsl_linalg_LU_decomp(matrix, perm, &s);
	gsl_linalg_LU_invert(matrix, perm, inverse);
	return inverse;
}

void matrixMultiply(gsl_matrix* A, gsl_matrix* B, gsl_matrix* C) {
	gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
                       1.0, A, B,
                       0.0, C);
	return;
}

void invertMatrix2x2(gsl_matrix* matrix, gsl_matrix* result) {
	double a = gsl_matrix_get(matrix, 0, 0);
	double b = gsl_matrix_get(matrix, 0, 1);
	double c = gsl_matrix_get(matrix, 1, 0);
	double d = gsl_matrix_get(matrix, 1, 1);
	double normalizer = 1.0f / (a*d - b*c);
	gsl_matrix_set(result, 0, 0, normalizer * d);
	gsl_matrix_set(result, 0, 1, normalizer * -b);
	gsl_matrix_set(result, 1, 0, normalizer * -c);
	gsl_matrix_set(result, 1, 1, normalizer * a);
	return;
}

void printMatrix(gsl_matrix* matrix) {
	int i;
	int j;
	int rows = matrix->size1;
	int cols = matrix->size2;
	for (i = 0; i < rows; i++) {
		for (j = 0; j < cols; j++) {
			printf("%2.4lf ", matrix->data[i * matrix->tda + j]);
		}
		printf("\n");
	}
	printf("\n");
}

void fillMatrix(gsl_matrix* matrix, double* data) {
	int i;
	int j;
	int rows = matrix->size1;
	int cols = matrix->size2;
	for (i = 0; i < rows; i++) {
		for (j = 0; j < cols; j++) {
			gsl_matrix_set(matrix, i, j, data[i * matrix->tda + j]);
		}
	}
	return;
}

